## Load packages
suppressMessages(library(tidyverse))
suppressMessages(library(lmerTest))
suppressMessages(library(MuMIn))
suppressMessages(library(broom.mixed))
suppressMessages(library(lavaan))
suppressMessages(library(reshape2))
suppressMessages(library(psych))
suppressMessages(library(glue))
suppressMessages(library(optimx))
# Read local files
source("utils.R")# Region reading time limits
rt.rel.sd = 2
rt.abs.min.pv = 300 # ms
rt.abs.max.pv = 3000 # ms
rt.abs.min.wm = 200 # ms
rt.abs.max.wm = 3000 # ms
rt.abs.max.irq = 120000 # ms
rt.abs.min.irq = 300 # ms
# Attention accuracy cutoff (retain if accuracy >= val)
cutoff.pv.accuracy = 0.66
cutoff.wm.accuracy = 0.66
cutoff.irq.attention = 0.66
# Exclude ppt with X% trials excluded
cutoff.removed.trials = 0.34
# lmer Control
optimxlmerControl <- lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))
optimxglmerControl <- glmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))# Load data
ppts.raw <- read.csv("data/vipv_participant.csv")
rating <- read.csv("data/vipv_rating.csv")
binary <- read.csv("data/vipv_binary.csv")
irq.og <- read.csv("data/irq.csv")# Preprocess participants
ppts.all <- preprocess.ppts(
ppts.raw, "pilot")
pv.all <- binary %>%
filter(participant_id %in% ppts.all$id,
task == "PV") %>%
mutate(accuracy = ifelse(is_correct == "True", 1, 0))
pv.critical.all <- pv.all %>% filter(item_type == "critical")
wm.all <- binary %>%
filter(participant_id %in% ppts.all$id,
task == "WM") %>%
mutate(accuracy = ifelse(is_correct == "True", 1, 0))
wm.critical.all <- wm.all %>% filter(item_type == "critical")
irq.all <- rating %>%
filter(participant_id %in% ppts.all$id,
scale == "IRQ")
rm(ppts.raw)
rm(binary)
rm(rating)# Exclude ppts
ppts.all <- exclude.ppts(ppts.all, native.eng=T, vision=T)ppts.all <- exclude.ppts.accuracy(ppts.all, pv.all, "pv", cutoff.pv.accuracy)95% of ppts got > 66% of PV trials correct.
ppts.all %>%
ggplot(aes(x = pv.accuracy, fill=excluded)) +
geom_histogram(bins=30) +
theme_minimal() +
scale_fill_manual("Excluded",
labels = c("Retained", "Excluded"),
values = c("#1eb809", "#cc0502")) +
labs(
x = "PV Accuracy",
y = "No. Participants"
)Participants who scored <66%:
ppts.all %>%
filter(pv.accuracy < 0.66) %>%
select(participant_id, pv.accuracy)| participant_id | pv.accuracy |
|---|---|
| 15 | 0.4739583 |
ppts.all <- exclude.ppts.accuracy(ppts.all, wm.all, "wm", cutoff.wm.accuracy)95% of ppts got > 66% of WM trials correct.
ppts.all %>%
ggplot(aes(x = wm.accuracy, fill=excluded)) +
scale_fill_manual("Excluded",
labels = c("Retained", "Excluded"),
values = c("#1eb809", "#cc0502")) +
geom_histogram(bins=30) +
theme_minimal() +
labs(
x = "WM Accuracy",
y = "No. Participants"
)Participants who scored <66%:
ppts.all %>%
filter(wm.accuracy < 0.66) %>%
select(participant_id, wm.accuracy)| participant_id | wm.accuracy |
|---|---|
| 15 | 0.4814815 |
ppts <- ppts.all %>% filter(excluded == 0)
summarise.exclusions(ppts.all)| Reason | Removed | (%) |
|---|---|---|
| ex.native_eng | 0 | 0.0 |
| ex.vision | 0 | 0.0 |
| ex.pv.accuracy | 1 | 4.8 |
| ex.wm.accuracy | 0 | 0.0 |
| —— | NA | NA |
| Total Removed | 1 | 4.8 |
| Retained | 20 | 95.2 |
Red dotted lines indicate exclusion criteria at 300ms and 3000ms.
pv.critical.all %>%
merge(ppts.all %>% select(participant_id, excluded)) %>%
ggplot(aes(x = reaction_time / 1000, fill=excluded)) +
geom_histogram(bins=30) +
scale_x_log10() +
theme_minimal() +
geom_vline(xintercept=rt.abs.min.pv / 1000, color="red", linetype="dashed") +
geom_vline(xintercept=rt.abs.max.pv / 1000, color="red", linetype="dashed") +
scale_fill_manual("Excluded",
labels = c("Retained", "Excluded"),
values = c("#1eb809", "#cc0502")) +
labs(
x = "Property verification reaction times (s)",
y = "No. Trials"
)pv.critical.all %>%
ggplot(aes(x = reaction_time, y=accuracy)) +
stat_summary_bin(fun.data = mean_cl_boot, geom="point") +
stat_summary_bin(fun.data = mean_cl_boot, geom="errorbar") +
geom_smooth(method="lm", formula="y~x") +
theme_minimal() +
geom_vline(xintercept=rt.abs.min.pv, color="red", linetype="dashed") +
geom_vline(xintercept=rt.abs.max.pv, color="red", linetype="dashed") +
scale_x_log10() +
labs()PV trial exclusion summary:
pv.critical.all <- pv.critical.all %>%
filter(participant_id %in% ppts$id, # Remove trials from excluded participants
) %>%
# Exclude extreme rts
exclude.relative.by.ppt(reaction_time, rt.rel.sd) %>%
exclude.absolute.by.ppt(reaction_time, rt.abs.min.pv, rt.abs.max.pv)
res = exclude.ppt.ex.trials(
ppts.all, pv.critical.all, cutoff.removed.trials, "ex.pv.ex.trials")
ppts.all <- res$ppts
pv.critical.all <- res$trials # TODO: Hacky. Better method?
pv.critical <- pv.critical.all %>% filter(excluded==FALSE)
summarise.exclusions(pv.critical.all)| Reason | Removed | (%) |
|---|---|---|
| ex.reaction_time.rel.lo | 0 | 0.0 |
| ex.reaction_time.rel.hi | 51 | 3.5 |
| ex.reaction_time.abs.lo | 0 | 0.0 |
| ex.reaction_time.abs.hi | 70 | 4.9 |
| ex.ppt.excluded | 0 | 0.0 |
| —— | NA | NA |
| Total Removed | 121 | 8.4 |
| Retained | 1319 | 91.6 |
Remove excluded participants:
ppts <- ppts.all %>% filter(excluded == 0)Red dotted lines indicate exclusion criteria at 200ms and 3000ms.
wm.critical.all %>%
merge(ppts.all %>% select(participant_id, excluded)) %>%
ggplot(aes(x = reaction_time / 1000, fill=excluded)) +
geom_histogram(bins=30) +
scale_x_log10() +
theme_minimal() +
geom_vline(xintercept=rt.abs.min.wm / 1000, color="red", linetype="dashed") +
geom_vline(xintercept=rt.abs.max.wm / 1000, color="red", linetype="dashed") +
scale_fill_manual("Excluded",
labels = c("Retained", "Excluded"),
values = c("#1eb809", "#cc0502")) +
labs(
x = "WM RT (s)",
y = "No. Trials"
)wm.critical.all %>%
ggplot(aes(x = reaction_time, y=accuracy)) +
stat_summary_bin(fun.data = mean_cl_boot, geom="point") +
stat_summary_bin(fun.data = mean_cl_boot, geom="errorbar") +
geom_smooth(method="lm", formula="y~x") +
theme_minimal() +
geom_vline(xintercept=rt.abs.min.wm, color="red", linetype="dashed") +
geom_vline(xintercept=rt.abs.max.wm, color="red", linetype="dashed") +
scale_x_log10() +
labs()WM trial exclusion summary:
wm.critical.all <- wm.critical.all %>%
filter(participant_id %in% ppts$id, # Remove trials from excluded participants
) %>%
# Exclude extreme rts
exclude.relative.by.ppt(reaction_time, rt.rel.sd) %>%
exclude.absolute.by.ppt(reaction_time, rt.abs.min.wm, rt.abs.max.wm)
res = exclude.ppt.ex.trials(
ppts.all, wm.critical.all, cutoff.removed.trials, "ex.wm.ex.trials")
ppts.all <- res$ppts
corsi.all <- res$trials
wm.critical <- wm.critical.all %>% filter(excluded==FALSE)
summarise.exclusions(wm.critical.all)| Reason | Removed | (%) |
|---|---|---|
| ex.reaction_time.rel.lo | 0 | 0.0 |
| ex.reaction_time.rel.hi | 50 | 5.2 |
| ex.reaction_time.abs.lo | 6 | 0.6 |
| ex.reaction_time.abs.hi | 11 | 1.1 |
| —— | NA | NA |
| Total Removed | 67 | 7.0 |
| Retained | 893 | 93.0 |
Remove excluded participants:
ppts <- ppts.all %>% filter(excluded == 0)pv.critical %>%
group_by(item) %>%
summarize(accuracy = mean(accuracy)) %>%
ggplot(aes(x = accuracy)) +
geom_histogram(aes(y=..count../sum(..count..) * 100), bins=10) +
labs(y = "Proportion of items")pv.critical %>%
group_by(item) %>%
summarize(accuracy = mean(accuracy)) %>%
filter(accuracy < 1) %>%
ggplot(aes(x = reorder(item, -accuracy), y = accuracy)) +
geom_point()# TODO: merge w/ stimuli data
pv.critical %>%
group_by(item) %>%
summarize(
reaction_time = mean(reaction_time),
accuracy = mean(accuracy)
) %>%
arrange(accuracy) %>%
head(10)| item | reaction_time | accuracy |
|---|---|---|
| 47 | 1626.000 | 0.1666667 |
| 65 | 1508.579 | 0.4736842 |
| 72 | 1873.353 | 0.5294118 |
| 69 | 1742.933 | 0.5333333 |
| 14 | 1506.222 | 0.5555556 |
| 17 | 1602.947 | 0.6315789 |
| 29 | 1749.947 | 0.6315789 |
| 70 | 1564.833 | 0.6666667 |
| 56 | 1847.250 | 0.6875000 |
| 36 | 1717.056 | 0.7222222 |
wm.dt <- wm.critical %>%
select(participant_id, trial_id, block_id, condition, version,
reaction_time, accuracy) %>%
rename(
wm.modality = condition,
wm.load = version,
wm.reaction_time = reaction_time,
wm.accuracy = accuracy
) %>%
mutate(
wm.load = factor(wm.load)
)
vipv <- merge(pv.critical %>% mutate(block_id = as.integer(block_id)), wm.dt, by=c("trial_id", "block_id", "participant_id"), all = F) %>%
rename (
pv.modality = condition,
pv.accuracy = accuracy,
pv.reaction_time = reaction_time,
) %>%
mutate(
pv.modality = factor(pv.modality, levels=c("vi", "au")),
wm.modality = factor(wm.modality, levels=c("vi", "au")),
acc.dt = pv.accuracy + wm.accuracy,
pv.rt.log = log(pv.reaction_time),
pv.rt.log.c = scale(pv.rt.log),
wm.rt.log = log(wm.reaction_time),
wm.rt.log.c = scale(wm.rt.log),
) %>%
arrange(trial_id)The raw relationship between RT and accuracy is negative, implying (to me) that there is no speed-accuracy trade-off(?). However, this could be due to a correlation between “easy” items (on which participants are fast and accurate) or “good” participants (who are fast and accurate).
vipv %>%
ggplot(aes(x = pv.reaction_time, y= pv.accuracy)) +
stat_summary_bin(geom="pointrange", bins=20, fun.data=mean_se) +
scale_x_log10() +
# geom_point() +
geom_smooth(method="lm", formula="y~x")## Warning: Removed 1 rows containing missing values (geom_segment).
We see the same trend within each condition combination
vipv %>%
ggplot(aes(x = pv.reaction_time, y= pv.accuracy, color=wm.load)) +
facet_grid(cols=vars(pv.modality), rows=vars(wm.modality), labeller="label_both") +
stat_summary_bin(geom="pointrange", bins=20, fun.data=mean_se) +
scale_x_log10() +
geom_smooth(method="lm", formula="y~x") +
theme_minimal()## Warning: Removed 4 rows containing missing values (geom_segment).
## Warning: Removed 5 rows containing missing values (geom_segment).
## Removed 5 rows containing missing values (geom_segment).
## Warning: Removed 2 rows containing missing values (geom_segment).
Participants who are faster are indeed more accurate
vipv %>%
group_by(participant_id) %>%
summarize(
pv.reaction_time = mean(pv.reaction_time),
pv.accuracy = mean(pv.accuracy)
) %>%
ggplot(aes(x = pv.reaction_time, y= pv.accuracy)) +
# stat_summary_bin(geom="pointrange", bins=20) +
scale_x_log10() +
geom_point() +
geom_smooth(method="lm", formula="y~x") +
labs(subtitle = "Grouped by participant")And accuracy and rt of items was also negatively correlated.
vipv %>%
group_by(item) %>%
summarize(
pv.reaction_time = mean(pv.reaction_time),
pv.accuracy = mean(pv.accuracy)
) %>%
ggplot(aes(x = pv.reaction_time, y= pv.accuracy)) +
# stat_summary_bin(geom="pointrange", bins=20) +
scale_x_log10() +
geom_point() +
geom_smooth(method="lm", formula="y~x") +
labs(subtitle = "Grouped by item")I tried to account for this by finding the by-ppt and by-item z-scores for accuracy and rt. Both continue to show a downward trend although potentially there’s a slight peak on the by-ppt plot.
vipv <- vipv %>%
group_by(item) %>%
mutate(
pv.accuracy.item.z = (pv.accuracy - mean(pv.accuracy)) / sd(pv.accuracy),
pv.rt.log.item.z = (pv.rt.log - mean(pv.rt.log)) / sd(pv.rt.log)
) %>%
ungroup() %>%
group_by(participant_id) %>%
mutate(
pv.accuracy.ppt.z = (pv.accuracy - mean(pv.accuracy)) / sd(pv.accuracy),
pv.rt.log.ppt.z = (pv.rt.log - mean(pv.rt.log)) / sd(pv.rt.log)
) %>%
ungroup() %>%
mutate(
pv.bis.ppt = pv.accuracy.ppt.z - pv.rt.log.ppt.z
)vipv %>%
ggplot(aes(x = pv.rt.log.ppt.z, y= pv.accuracy.ppt.z)) +
stat_summary_bin(geom="pointrange", bins=20, fun.data=mean_se) +
# geom_point() +
geom_smooth(method="lm", formula="y~x")## Warning: Removed 62 rows containing non-finite values (stat_summary_bin).
## Warning: Removed 62 rows containing non-finite values (stat_smooth).
## Warning: Removed 2 rows containing missing values (geom_segment).
vipv %>%
ungroup() %>%
ggplot(aes(x = pv.rt.log.item.z, y= pv.accuracy.item.z)) +
stat_summary_bin(geom="pointrange", bins=20) +
# geom_point() +
geom_smooth(method="lm", formula="y~x")## Warning: Removed 476 rows containing non-finite values (stat_summary_bin).
## No summary function supplied, defaulting to `mean_se()`
## Warning: Removed 476 rows containing non-finite values (stat_smooth).
## Warning: Removed 2 rows containing missing values (geom_segment).
For PV RT we see the opposite cross-over interaction in the high-load condition, compared to the one observed by Vermeulen et al. Ppts are faster to respond to trials where the modality of the WM and PV trials match.
vipv %>%
filter(wm.accuracy == 1,
pv.accuracy == 1) %>%
ggplot(aes(x=pv.modality, y = pv.rt.log, color=as.factor(wm.modality))) +
stat_summary(fun="mean", geom="point",size=5, position=position_dodge(width=1)) +
stat_summary(fun.data="mean_se", geom="errorbar", position=position_dodge(width=1)) +
facet_grid(cols=vars(wm.load), labeller = label_both) +
theme(
) +
labs(
x = "PV Modality",
y = "PV Reaction Time (log)",
color = "WM Modality"
)The 3-way interaction has a non-sigificant, negative effect on RT’s.
m.vipv.pv.rt.base <- lmer(
pv.rt.log ~ pv.modality * wm.modality +
pv.modality * wm.load + wm.modality * wm.load +
(1 | participant_id) + (1 | item_id),
data = vipv %>%
filter(wm.accuracy == 1,
pv.accuracy == 1), REML=F)
m.vipv.pv.rt <- lmer(
pv.rt.log ~ pv.modality * wm.modality * wm.load +
(1 | participant_id) + (1 | item_id),
data = vipv %>%
filter(wm.accuracy == 1,
pv.accuracy == 1), REML=F)
summary(m.vipv.pv.rt)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula:
## pv.rt.log ~ pv.modality * wm.modality * wm.load + (1 | participant_id) +
## (1 | item_id)
## Data: vipv %>% filter(wm.accuracy == 1, pv.accuracy == 1)
##
## AIC BIC logLik deviance df.resid
## 96.4 146.5 -37.2 74.4 689
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1724 -0.6633 -0.0587 0.5589 3.1892
##
## Random effects:
## Groups Name Variance Std.Dev.
## item_id (Intercept) 0.01384 0.1176
## participant_id (Intercept) 0.03180 0.1783
## Residual 0.05383 0.2320
## Number of obs: 700, groups: item_id, 48; participant_id, 20
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 7.222696 0.053352 53.139006 135.379
## pv.modalityau 0.021376 0.049588 120.453682 0.431
## wm.modalityau 0.030580 0.036353 651.247657 0.841
## wm.load3 -0.041210 0.036728 653.536613 -1.122
## pv.modalityau:wm.modalityau -0.047143 0.051382 652.189901 -0.918
## pv.modalityau:wm.load3 0.008185 0.051268 650.476869 0.160
## wm.modalityau:wm.load3 0.031695 0.051444 658.292137 0.616
## pv.modalityau:wm.modalityau:wm.load3 -0.041058 0.073240 656.550607 -0.561
## Pr(>|t|)
## (Intercept) <2e-16 ***
## pv.modalityau 0.667
## wm.modalityau 0.401
## wm.load3 0.262
## pv.modalityau:wm.modalityau 0.359
## pv.modalityau:wm.load3 0.873
## wm.modalityau:wm.load3 0.538
## pv.modalityau:wm.modalityau:wm.load3 0.575
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pv.mdl wm.mdl wm.ld3 pv.m:. pv.:.3 wm.:.3
## pv.modality -0.475
## wm.modality -0.355 0.383
## wm.load3 -0.347 0.375 0.518
## pv.mdlty:w. 0.252 -0.516 -0.709 -0.370
## pv.mdlty:.3 0.250 -0.508 -0.373 -0.719 0.504
## wm.mdlty:.3 0.253 -0.274 -0.710 -0.732 0.505 0.526
## pv.mdl:.:.3 -0.179 0.362 0.501 0.518 -0.706 -0.720 -0.705
a.vipv.pv.rt <- anova(m.vipv.pv.rt.base, m.vipv.pv.rt)
a.vipv.pv.rt| npar | AIC | BIC | logLik | deviance | Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|---|
| m.vipv.pv.rt.base | 10 | 94.70508 | 140.2159 | -37.35254 | 74.70508 | NA | NA | NA |
| m.vipv.pv.rt | 11 | 96.39114 | 146.4530 | -37.19557 | 74.39114 | 0.3139382 | 1 | 0.5752735 |
PV accuracy appears to show a 2-way interaction in the expected direction, with accuracy generally lower for same-modality trials in both the low and high-load conditions.
vipv %>%
filter(wm.accuracy == 1) %>%
ggplot(aes(x=pv.modality, y = pv.accuracy, color=as.factor(wm.modality))) +
stat_summary(fun="mean", geom="point",size=5, position=position_dodge(width=1)) +
stat_summary(fun.data="mean_se", geom="errorbar", position=position_dodge(width=1)) +
facet_grid(cols=vars(wm.load), labeller = label_value) +
theme(
) +
labs(
x = "PV Modality",
y = "PV Accuracy",
color = "WM Modality"
)A linear model shows no 2-way interaction (or any other effects).
m.vipv.pv.acc.base <- glmer(
pv.accuracy ~ 1 + pv.modality + wm.modality + wm.load +
pv.modality:wm.modality + wm.modality:wm.load + pv.modality:wm.load +
(1 | participant_id) + (1 | item_id),
data = vipv %>%
filter(wm.accuracy == 1),
family="binomial", control=optimxglmerControl)
m.vipv.pv.acc <- glmer(
pv.accuracy ~ 1 + pv.modality + wm.modality + wm.load +
pv.modality:wm.modality + wm.modality:wm.load + pv.modality:wm.load +
pv.modality:wm.modality:wm.load +
(1 | participant_id) + (1 | item_id),
data = vipv %>%
filter(wm.accuracy == 1),
family="binomial", control=optimxglmerControl)
summary(m.vipv.pv.acc)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## pv.accuracy ~ 1 + pv.modality + wm.modality + wm.load + pv.modality:wm.modality +
## wm.modality:wm.load + pv.modality:wm.load + pv.modality:wm.modality:wm.load +
## (1 | participant_id) + (1 | item_id)
## Data: vipv %>% filter(wm.accuracy == 1)
## Control: optimxglmerControl
##
## AIC BIC logLik deviance df.resid
## 410.4 457.0 -195.2 390.4 776
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.8547 0.0450 0.0975 0.1880 3.2077
##
## Random effects:
## Groups Name Variance Std.Dev.
## item_id (Intercept) 7.882 2.808
## participant_id (Intercept) 1.426 1.194
## Number of obs: 786, groups: item_id, 48; participant_id, 20
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.70721 0.98846 4.762 1.92e-06 ***
## pv.modalityau -0.76484 1.11833 -0.684 0.494
## wm.modalityau -0.02312 0.66705 -0.035 0.972
## wm.load3 -0.27648 0.63688 -0.434 0.664
## pv.modalityau:wm.modalityau 0.13207 0.89013 0.148 0.882
## wm.modalityau:wm.load3 0.96221 0.98268 0.979 0.327
## pv.modalityau:wm.load3 0.86039 0.88375 0.974 0.330
## pv.modalityau:wm.modalityau:wm.load3 -1.09267 1.31435 -0.831 0.406
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pv.mdl wm.mdl wm.ld3 pv.m:. wm.:.3 pv.:.3
## pv.modality -0.589
## wm.modality -0.339 0.291
## wm.load3 -0.352 0.299 0.488
## pv.mdlty:w. 0.266 -0.389 -0.757 -0.367
## wm.mdlty:.3 0.281 -0.213 -0.695 -0.682 0.517
## pv.mdlty:.3 0.286 -0.380 -0.363 -0.729 0.500 0.505
## pv.mdl:.:.3 -0.224 0.270 0.531 0.518 -0.675 -0.760 -0.704
a.vipv.pv.acc <- anova(m.vipv.pv.acc.base, m.vipv.pv.acc)
a.vipv.pv.acc| npar | AIC | BIC | logLik | deviance | Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|---|
| m.vipv.pv.acc.base | 9 | 409.0486 | 451.0512 | -195.5243 | 391.0486 | NA | NA | NA |
| m.vipv.pv.acc | 10 | 410.3519 | 457.0215 | -195.1760 | 390.3519 | 0.6967111 | 1 | 0.403891 |
In an attempt to combine RT and Accuracy, I used the Balanced Integration of Scores (BIS) (Liesefeld & Janczyk, 2018), essentially ppt-wise z(acc) - z(rt).
It seems to mostly be driven by RT, showing the same rough pattern of results (better performance on same-modality trials).
vipv %>%
filter(wm.accuracy == 1) %>%
ggplot(aes(x=pv.modality, y = pv.bis.ppt, color=as.factor(wm.modality))) +
stat_summary(fun="mean", geom="point",size=5, position=position_dodge(width=1)) +
stat_summary(fun.data="mean_se", geom="errorbar", position=position_dodge(width=1)) +
facet_grid(cols=vars(wm.load), labeller = label_both) +
theme(
) +
labs(
x = "PV Modality",
y = "PV BIS",
color = "WM Modality"
)## Warning: Removed 57 rows containing non-finite values (stat_summary).
## Removed 57 rows containing non-finite values (stat_summary).
No significant effects.
m.vipv.pv.bis.base <- lmer(
pv.bis.ppt ~ 1 + pv.modality + wm.modality + wm.load +
pv.modality:wm.modality + wm.modality:wm.load + pv.modality:wm.load +
(1 | participant_id) + (1 | item_id),
data = vipv %>%
filter(wm.accuracy == 1), REML=F)## boundary (singular) fit: see help('isSingular')
m.vipv.pv.bis <- lmer(
pv.bis.ppt ~ 1 + pv.modality + wm.modality + wm.load +
pv.modality:wm.modality + wm.modality:wm.load + pv.modality:wm.load +
pv.modality:wm.modality:wm.load +
(1 | participant_id) + (1 | item_id),
data = vipv %>%
filter(wm.accuracy == 1), REML=F)## boundary (singular) fit: see help('isSingular')
summary(m.vipv.pv.bis)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula:
## pv.bis.ppt ~ 1 + pv.modality + wm.modality + wm.load + pv.modality:wm.modality +
## wm.modality:wm.load + pv.modality:wm.load + pv.modality:wm.modality:wm.load +
## (1 | participant_id) + (1 | item_id)
## Data: vipv %>% filter(wm.accuracy == 1)
##
## AIC BIC logLik deviance df.resid
## 2476.1 2526.6 -1227.0 2454.1 718
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1167 -0.5397 0.1048 0.6036 3.6015
##
## Random effects:
## Groups Name Variance Std.Dev.
## item_id (Intercept) 0.8319 0.9121
## participant_id (Intercept) 0.0000 0.0000
## Residual 1.4617 1.2090
## Number of obs: 729, groups: item_id, 48; participant_id, 18
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 0.074414 0.228992 86.665046 0.325
## pv.modalityau -0.285895 0.319275 82.041512 -0.895
## wm.modalityau -0.067562 0.187214 688.335785 -0.361
## wm.load3 0.171609 0.186248 688.128275 0.921
## pv.modalityau:wm.modalityau 0.007266 0.258522 688.884962 0.028
## wm.modalityau:wm.load3 -0.030716 0.263751 691.066143 -0.116
## pv.modalityau:wm.load3 0.034346 0.258473 687.458494 0.133
## pv.modalityau:wm.modalityau:wm.load3 0.134498 0.368936 689.849370 0.365
## Pr(>|t|)
## (Intercept) 0.746
## pv.modalityau 0.373
## wm.modalityau 0.718
## wm.load3 0.357
## pv.modalityau:wm.modalityau 0.978
## wm.modalityau:wm.load3 0.907
## pv.modalityau:wm.load3 0.894
## pv.modalityau:wm.modalityau:wm.load3 0.716
##
## Correlation of Fixed Effects:
## (Intr) pv.mdl wm.mdl wm.ld3 pv.m:. wm.:.3 pv.:.3
## pv.modality -0.717
## wm.modality -0.421 0.302
## wm.load3 -0.417 0.299 0.518
## pv.mdlty:w. 0.305 -0.403 -0.724 -0.375
## wm.mdlty:.3 0.300 -0.215 -0.713 -0.723 0.516
## pv.mdlty:.3 0.300 -0.397 -0.373 -0.721 0.502 0.521
## pv.mdl:.:.3 -0.215 0.282 0.509 0.517 -0.702 -0.715 -0.715
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
a.vipv.pv.bis <- anova(m.vipv.pv.bis.base, m.vipv.pv.bis)
a.vipv.pv.bis| npar | AIC | BIC | logLik | deviance | Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|---|
| m.vipv.pv.bis.base | 10 | 2474.202 | 2520.119 | -1227.101 | 2454.202 | NA | NA | NA |
| m.vipv.pv.bis | 11 | 2476.070 | 2526.578 | -1227.035 | 2454.070 | 0.1328024 | 1 | 0.7155439 |
WM RTs also show what could be a crossover in the high-load condition, again with faster RTs for same-modality trials.
vipv %>%
filter(wm.accuracy == 1,
pv.accuracy == 1) %>%
ggplot(aes(x=pv.modality, y = wm.rt.log, color=as.factor(wm.modality))) +
stat_summary(fun="mean", geom="point",size=5, position=position_dodge(width=1)) +
stat_summary(fun.data="mean_se", geom="errorbar", position=position_dodge(width=1)) +
facet_grid(cols=vars(wm.load), labeller = label_value) +
theme(
) +
labs(
x = "PV Modality",
y = "WM Reaction Time (log s)",
color = "WM Modality"
)The 3-way interaction has a negative effect which is not significant (p=0.23).
m.vipv.wm.rt.base <- lmer(
wm.rt.log ~ 1 + pv.modality + wm.modality + wm.load +
pv.modality:wm.modality + wm.modality:wm.load + pv.modality:wm.load +
(1 | participant_id) + (1 | item_id),
data = vipv %>%
filter(wm.accuracy == 1,
pv.accuracy == 1), REML=F)## boundary (singular) fit: see help('isSingular')
m.vipv.wm.rt <- lmer(
wm.rt.log ~ 1 + pv.modality + wm.modality + wm.load +
pv.modality:wm.modality + wm.modality:wm.load + pv.modality:wm.load +
pv.modality:wm.modality:wm.load +
(1 | participant_id) + (1 | item_id),
data = vipv %>%
filter(wm.accuracy == 1,
pv.accuracy == 1), REML=F)## boundary (singular) fit: see help('isSingular')
summary(m.vipv.wm.rt)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula:
## wm.rt.log ~ 1 + pv.modality + wm.modality + wm.load + pv.modality:wm.modality +
## wm.modality:wm.load + pv.modality:wm.load + pv.modality:wm.modality:wm.load +
## (1 | participant_id) + (1 | item_id)
## Data: vipv %>% filter(wm.accuracy == 1, pv.accuracy == 1)
##
## AIC BIC logLik deviance df.resid
## 860.6 910.7 -419.3 838.6 689
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3806 -0.6801 -0.0311 0.6325 3.2576
##
## Random effects:
## Groups Name Variance Std.Dev.
## item_id (Intercept) 0.00000 0.0000
## participant_id (Intercept) 0.03897 0.1974
## Residual 0.18256 0.4273
## Number of obs: 700, groups: item_id, 48; participant_id, 20
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 6.536556 0.064548 66.203990 101.267
## pv.modalityau 0.014769 0.064593 680.516201 0.229
## wm.modalityau -0.032284 0.064999 680.021736 -0.497
## wm.load3 -0.060652 0.065646 680.663989 -0.924
## pv.modalityau:wm.modalityau -0.003506 0.091743 680.680192 -0.038
## wm.modalityau:wm.load3 0.048504 0.091003 679.973265 0.533
## pv.modalityau:wm.load3 0.050945 0.091976 681.357950 0.554
## pv.modalityau:wm.modalityau:wm.load3 -0.158006 0.130061 681.253933 -1.215
## Pr(>|t|)
## (Intercept) <2e-16 ***
## pv.modalityau 0.819
## wm.modalityau 0.620
## wm.load3 0.356
## pv.modalityau:wm.modalityau 0.970
## wm.modalityau:wm.load3 0.594
## pv.modalityau:wm.load3 0.580
## pv.modalityau:wm.modalityau:wm.load3 0.225
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pv.mdl wm.mdl wm.ld3 pv.m:. wm.:.3 pv.:.3
## pv.modality -0.531
## wm.modality -0.528 0.527
## wm.load3 -0.523 0.525 0.519
## pv.mdlty:w. 0.375 -0.706 -0.709 -0.371
## wm.mdlty:.3 0.377 -0.378 -0.714 -0.721 0.508
## pv.mdlty:.3 0.375 -0.705 -0.371 -0.717 0.499 0.516
## pv.mdl:.:.3 -0.266 0.500 0.501 0.508 -0.707 -0.702 -0.709
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
a.vipv.wm.rt <- anova(m.vipv.wm.rt.base, m.vipv.wm.rt)
a.vipv.wm.rt| npar | AIC | BIC | logLik | deviance | Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|---|
| m.vipv.wm.rt.base | 10 | 860.0706 | 905.5814 | -420.0353 | 840.0706 | NA | NA | NA |
| m.vipv.wm.rt | 11 | 860.5966 | 910.6585 | -419.2983 | 838.5966 | 1.474016 | 1 | 0.2247131 |
WM Accuracy seems to show a main effect of load (which makes sense).
vipv %>%
filter(pv.accuracy == 1) %>%
ggplot(aes(x=pv.modality, y = wm.accuracy, color=as.factor(wm.modality))) +
stat_summary(fun="mean", geom="point",size=5, position=position_dodge(width=1)) +
stat_summary(fun.data="mean_se", geom="errorbar", position=position_dodge(width=1)) +
facet_grid(cols=vars(wm.load), labeller = label_value) +
theme(
) +
labs(
x = "PV Modality",
y = "WM Accuracy",
color = "WM Modality"
)The load effect, however, is nonsignificant (p=0.36)
m.vipv.wm.acc.base <- glmer(
wm.accuracy ~ 1 + pv.modality + wm.modality + wm.load +
pv.modality:wm.modality + wm.modality:wm.load + pv.modality:wm.load +
(1 | participant_id) + (1 | item_id),
data = vipv %>%
filter(pv.accuracy == 1),
family="binomial", control=optimxglmerControl)
m.vipv.wm.acc <- glmer(
wm.accuracy ~ 1 + pv.modality + wm.modality + wm.load +
pv.modality:wm.modality + wm.modality:wm.load + pv.modality:wm.load +
pv.modality:wm.modality:wm.load +
(1 | participant_id) + (1 | item_id),
data = vipv %>%
filter(pv.accuracy == 1),
family="binomial", control=optimxglmerControl)
summary(m.vipv.wm.acc)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## wm.accuracy ~ 1 + pv.modality + wm.modality + wm.load + pv.modality:wm.modality +
## wm.modality:wm.load + pv.modality:wm.load + pv.modality:wm.modality:wm.load +
## (1 | participant_id) + (1 | item_id)
## Data: vipv %>% filter(pv.accuracy == 1)
## Control: optimxglmerControl
##
## AIC BIC logLik deviance df.resid
## 291.7 337.7 -135.9 271.7 727
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.2790 0.1175 0.1613 0.2319 0.5676
##
## Random effects:
## Groups Name Variance Std.Dev.
## item_id (Intercept) 0.5681 0.7537
## participant_id (Intercept) 0.5692 0.7544
## Number of obs: 737, groups: item_id, 48; participant_id, 20
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.75798 0.68542 5.483 4.19e-08 ***
## pv.modalityau 0.35789 0.90644 0.395 0.693
## wm.modalityau 0.59758 0.95432 0.626 0.531
## wm.load3 -0.68847 0.75335 -0.914 0.361
## pv.modalityau:wm.modalityau 0.10432 1.53791 0.068 0.946
## wm.modalityau:wm.load3 0.03283 1.18269 0.028 0.978
## pv.modalityau:wm.load3 -0.97142 1.05270 -0.923 0.356
## pv.modalityau:wm.modalityau:wm.load3 -0.36073 1.78052 -0.203 0.839
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pv.mdl wm.mdl wm.ld3 pv.m:. wm.:.3 pv.:.3
## pv.modality -0.589
## wm.modality -0.543 0.435
## wm.load3 -0.721 0.533 0.513
## pv.mdlty:w. 0.313 -0.558 -0.627 -0.317
## wm.mdlty:.3 0.460 -0.343 -0.803 -0.646 0.498
## pv.mdlty:.3 0.460 -0.811 -0.381 -0.711 0.487 0.461
## pv.mdl:.:.3 -0.271 0.488 0.543 0.427 -0.866 -0.665 -0.604
a.vipv.wm.acc <- anova(m.vipv.wm.acc.base, m.vipv.wm.acc)
a.vipv.wm.acc| npar | AIC | BIC | logLik | deviance | Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|---|
| m.vipv.wm.acc.base | 9 | 289.7638 | 331.1871 | -135.8819 | 271.7638 | NA | NA | NA |
| m.vipv.wm.acc | 10 | 291.7225 | 337.7483 | -135.8612 | 271.7225 | 0.0413126 | 1 | 0.8389357 |
p.val <- c(a.vipv.pv.rt$`Pr(>Chisq)`[2],
a.vipv.pv.acc$`Pr(>Chisq)`[2],
a.vipv.wm.rt$`Pr(>Chisq)`[2],
a.vipv.wm.acc$`Pr(>Chisq)`[2])
test <- c("PV RT", "PV Acc", "WM RT", "WM Acc")
p.val.adj <- stats::p.adjust(p.val, method="BH")
data.frame(test, p.val, p.val.adj)| test | p.val | p.val.adj |
|---|---|---|
| PV RT | 0.5752735 | 0.7670314 |
| PV Acc | 0.4038910 | 0.7670314 |
| WM RT | 0.2247131 | 0.7670314 |
| WM Acc | 0.8389357 | 0.8389357 |
irq.attention.answers <- data.frame(
item_id = c(37, 38, 39),
correct_response = c(1,5,1))
irq.all <- irq.all %>%
mutate(item_type = recode(item_type, catch = "attention"))
ppts.all <- exclude.ppts.attention(ppts.all, irq.all, irq.attention.answers, "irq", cutoff.irq.attention)100% of ppts got > 66% IRQ attention checks correct.
ppts.all %>%
ggplot(aes(x = irq.attention.accuracy)) +
geom_histogram() +
theme_minimal() +
labs(
x = "IRQ Attention Accuracy",
y = "No. Participants"
)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Participants who failed IRQ attention check.
ppts.all %>%
filter(irq.attention.accuracy < 0.66) %>%
select(participant_id, irq.attention.accuracy)| participant_id | irq.attention.accuracy |
|---|
IRQ Accuracy by question
irq.all %>%
filter(item_type == "attention") %>%
merge(irq.attention.answers) %>%
mutate(accuracy = ifelse(response == correct_response, 1, 0)) %>%
group_by(item_id) %>%
summarize(accuracy = mean(accuracy), .groups="drop")| item_id | accuracy |
|---|---|
| 37 | 1 |
| 38 | 1 |
| 39 | 1 |
irq.all <- irq.all %>%
filter(participant_id %in% ppts$id, # Remove trials from excluded participants
item_type != "attention" # Remove catch trials
) %>%
# Exclude extreme rts
# exclude.relative.by.ppt(reaction_time, rt.rel.sd) %>%
exclude.absolute.by.ppt(reaction_time, rt.abs.min.irq, rt.abs.max.irq)
res = exclude.ppt.ex.trials(
ppts.all, irq.all, cutoff.removed.trials, "ex.irq.ex.trials")
ppts.all <- res$ppts
irq.all <- res$trials
irq <- irq.all %>% filter(excluded==FALSE)
summarise.exclusions(irq.all)| Reason | Removed | (%) |
|---|---|---|
| ex.reaction_time.abs.lo | 0 | 0 |
| ex.reaction_time.abs.hi | 0 | 0 |
| ex.ppt.excluded | 0 | 0 |
| —— | NA | NA |
| Total Removed | 0 | 0 |
| Retained | 720 | 100 |
Pivot the IRQ question data into (item x ppt response)
irq <- irq.all %>%
filter(participant_id %in% ppts$id)
irq.qs <- irq %>%
select(item_id, response, participant_id) %>%
pivot_wider(names_from=item_id, values_from=response) %>%
arrange(participant_id) %>%
select(-participant_id) %>%
select(order(as.numeric(colnames(.))))
# data.frame()
# Reverse loadings for reversed q's
irq.qs[10] = 6 - irq.qs[10]
irq.qs[19] = 6 - irq.qs[19]
irq.qs[33] = 6 - irq.qs[33]
# irq.qsParallel analysis recommends 2 factors.
fa.parallel(irq.qs)## In smc, smcs < 0 were set to .0
## In smc, smcs < 0 were set to .0
## In smc, smcs < 0 were set to .0
## Warning in cor.smooth(r): Matrix was not positive definite, smoothing was done
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.
## In factor.scores, the correlation matrix is singular, the pseudo inverse is used
## I was unable to calculate the factor score weights, factor loadings used instead
## Parallel analysis suggests that the number of factors = 2 and the number of components = 2
Re-performing FA with 4 factors seems to produce similar-looking loadings to the original paper (in the same order, i.e. MR1=verbal, MR2=orthographics, MR3=visual, MR4=manipulation)
irq.fa.4 <- psych::fa(irq.qs, nfactors=4,
rotate="oblimin", # Oblique rotation as in original paper
fm="minres" # Default minimum residual factoring method
)## In smc, smcs < 0 were set to .0
## In smc, smcs < 0 were set to .0
## In smc, smcs < 0 were set to .0
## Loading required namespace: GPArotation
## Warning in cor.smooth(r): Matrix was not positive definite, smoothing was done
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.
## In factor.scores, the correlation matrix is singular, the pseudo inverse is used
## I was unable to calculate the factor score weights, factor loadings used instead
loadings <- irq.fa.4$loadings
loadings <- as.data.frame(loadings[1:36,1:4])
loadings$item_id <- row.names(loadings)
loadings$item.factor <- c(rep("Visual", 10), rep("Verbal", 12), rep("Orthographic", 6), rep("Manipulation", 8))
loadings.pivot <- loadings %>%
pivot_longer(cols=c(MR1, MR2, MR3, MR4),
names_to="Factor", values_to="Loading New")
loadings.pivot %>%
ggplot(aes(x = Factor, y = `Loading New`, color=item.factor)) +
stat_summary(fun.data=mean_se, geom="pointrange") # geom_point(stat="identity", position="dodge")
# loadingsCompare loadings for each item between original and new data analysis. The match looks fairly good visually, the same set of items in the high end of each axis.
colnames(loadings) <- c("Verbal", "Orthographic", "Visual", "Manipulation", "item_id")
irq.og$item.factor <- c(rep("Visual", 10), rep("Verbal", 12), rep("Orthographic", 6), rep("Manipulation", 8), rep("Filler", 3))
irq.og.loadings <- irq.og %>%
filter(item_type == "critical") %>% # remove fillers
select(item_id, item.factor, Visual, Verbal, Orthographic, Manipulation)
loadings.pivot <- loadings %>%
select(Visual, Verbal, Orthographic, Manipulation, item_id) %>%
pivot_longer(cols=c(Visual, Verbal, Orthographic, Manipulation),
names_to="Factor", values_to="Loading New")
loadings.og.pivot <- irq.og.loadings %>%
pivot_longer(cols=c(Visual, Verbal, Orthographic, Manipulation),
names_to="Factor", values_to="Loading OG")
loadings.all <- merge(loadings.pivot, loadings.og.pivot)
loadings.all %>%
ggplot(aes(x = `Loading OG`, y = `Loading New`, color=item.factor)) +
geom_point() +
geom_text(aes(label=item_id), nudge_x = 0.03, nudge_y = 0.03) +
facet_wrap(facets=vars(Factor)) +
theme_minimal()Merge IRQ data onto ppts:
I’m using EFA fit as it seems to be what was used in the original paper.
# Merge ppt scores
ppt.scores <- as.data.frame(irq.fa.4$scores)
colnames(ppt.scores) <- c("Verbal", "Orthographic", "Visual", "Manipulation")
ppts <- cbind(ppts,ppt.scores)m.vipv.pv.acc.ixn <- glmer(
pv.accuracy ~ 1 + pv.modality + wm.modality + wm.load +
pv.modality:wm.modality + wm.modality:wm.load + pv.modality:wm.load +
pv.modality:wm.modality:wm.load +
(
pv.modality + wm.modality + wm.load +
pv.modality:wm.modality + wm.modality:wm.load + pv.modality:wm.load +
pv.modality:wm.modality:wm.load | participant_id) + (1 | item_id),
data = vipv %>% filter(wm.accuracy == 1),
family="binomial", control=optimxglmerControl)## Warning in optwrap(optimizer, devfun, start, rho$lower, control = control, :
## convergence code 1 from optimx: none
## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## Warning in optwrap(optimizer, devfun, start, rho$lower, control = control, :
## convergence code 1 from optimx: none
## boundary (singular) fit: see help('isSingular')
summary(m.vipv.pv.acc.ixn)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## pv.accuracy ~ 1 + pv.modality + wm.modality + wm.load + pv.modality:wm.modality +
## wm.modality:wm.load + pv.modality:wm.load + pv.modality:wm.modality:wm.load +
## (pv.modality + wm.modality + wm.load + pv.modality:wm.modality +
## wm.modality:wm.load + pv.modality:wm.load + pv.modality:wm.modality:wm.load |
## participant_id) + (1 | item_id)
## Data: vipv %>% filter(wm.accuracy == 1)
## Control: optimxglmerControl
##
## AIC BIC logLik deviance df.resid
## 466.1 676.1 -188.0 376.1 741
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5294 0.0259 0.0676 0.1542 2.1911
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## item_id (Intercept) 10.940 3.308
## participant_id (Intercept) 6.002 2.450
## pv.modalityau 4.211 2.052 -0.87
## wm.modalityau 7.817 2.796 -0.77
## wm.load3 8.404 2.899 -0.92
## pv.modalityau:wm.modalityau 4.032 2.008 0.66
## wm.modalityau:wm.load3 16.235 4.029 0.85
## pv.modalityau:wm.load3 14.810 3.848 0.96
## pv.modalityau:wm.modalityau:wm.load3 23.197 4.816 -0.57
##
##
##
##
## 0.95
## 0.87 0.82
## -0.88 -0.80 -0.54
## -0.98 -0.93 -0.78 0.94
## -0.94 -0.81 -0.93 0.76 0.90
## 0.72 0.56 0.36 -0.94 -0.82 -0.67
## Number of obs: 786, groups: item_id, 48; participant_id, 20
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.1635 1.6807 3.667 0.000245 ***
## pv.modalityau -1.7539 1.7290 -1.014 0.310391
## wm.modalityau -0.4129 1.4927 -0.277 0.782046
## wm.load3 -1.2711 1.3769 -0.923 0.355930
## pv.modalityau:wm.modalityau 0.9880 1.6144 0.612 0.540524
## wm.modalityau:wm.load3 2.6801 2.3338 1.148 0.250812
## pv.modalityau:wm.load3 2.6345 1.8241 1.444 0.148663
## pv.modalityau:wm.modalityau:wm.load3 -3.5109 2.7422 -1.280 0.200427
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pv.mdl wm.mdl wm.ld3 pv.m:. wm.:.3 pv.:.3
## pv.modality -0.719
## wm.modality -0.569 0.531
## wm.load3 -0.674 0.603 0.676
## pv.mdlty:w. 0.511 -0.556 -0.860 -0.563
## wm.mdlty:.3 0.432 -0.393 -0.675 -0.585 0.598
## pv.mdlty:.3 0.591 -0.581 -0.565 -0.829 0.574 0.530
## pv.mdl:.:.3 -0.363 0.373 0.546 0.461 -0.639 -0.867 -0.626
## optimizer (optimx) convergence code: 1 (none)
## boundary (singular) fit: see help('isSingular')
## Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
We take the random slope for the domainVisual effect for each participant in a linear model.
ppt.ranef = ranef(m.vipv.pv.acc.ixn)$participant_id
ppt.ranef <- ppt.ranef %>%
mutate(participant_id = rownames(ppt.ranef))
ppt.ranef$ixn.full = ppt.ranef$`pv.modalityau:wm.modalityau:wm.load3` + fixef(m.vipv.pv.acc.ixn)["pv.modalityau:wm.modalityau:wm.load3"]There is a some variance in the effect by participant, with participants varying from -6.1653481 to 4.1852347.
ppt.ranef %>%
ggplot(aes(x = reorder(participant_id, -ixn.full), y = ixn.full)) +
geom_hline(yintercept=fixef(m.vipv.pv.acc.ixn)["pv.modalityau:wm.modalityau:wm.load3"], color="red", linetype="dashed") +
geom_point() +
theme_minimal() +
theme(
axis.text.x = element_blank()
) +
labs(
subtitle = "Participant Random Slopes for 3-way Interaction Effect",
title = "PV Accuracy Random Slopes",
y = "3-way Interaction Effect (Fixed Effect + Random Slope)",
x = "Participant Id",
caption = "Red dashed line represents fixed effect"
)ppts <- merge(ppts, ppt.ranef, by="participant_id", all.y=F)Correlations of IRQ w/ ppt random-slopes don’t show much of interest.
ppts %>%
pivot_longer(cols=c(Visual, Verbal, Orthographic, Manipulation), names_to="factor", values_to="score") %>%
ggplot(aes(x = score, y = ixn.full)) +
geom_point() +
geom_smooth(method="lm", formula="y~x") +
facet_wrap(facets=vars(factor)) +
theme_minimal() +
labs(
# title="IRQ vs PR Domain ∆",
# subtitle = "By-Participant",
# x = "IRQ Dimension Score",
# y = "PR Accuracy Domain ∆ (Visual - Social)"
)Simple linear models show no effect of either variable.
summary(lm(ixn.full ~ Visual, data=ppts))
summary(lm(ixn.full ~ Manipulation, data=ppts))##
## Call:
## lm(formula = ixn.full ~ Visual, data = ppts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3426 -1.8269 -0.6854 0.6660 6.8438
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.91471 0.66079 -4.411 0.000337 ***
## Visual 0.06123 0.17879 0.342 0.735958
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.955 on 18 degrees of freedom
## Multiple R-squared: 0.006474, Adjusted R-squared: -0.04872
## F-statistic: 0.1173 on 1 and 18 DF, p-value: 0.736
##
##
## Call:
## lm(formula = ixn.full ~ Manipulation, data = ppts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3340 -1.9139 -0.6915 0.5977 7.1268
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.91471 0.66278 -4.398 0.000347 ***
## Manipulation -0.01703 0.18418 -0.092 0.927362
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.964 on 18 degrees of freedom
## Multiple R-squared: 0.0004746, Adjusted R-squared: -0.05505
## F-statistic: 0.008547 on 1 and 18 DF, p-value: 0.9274
We hypothesized that verbalizers would show a bigger main effect of WM modality, because auditory info would interfere with their default representation strategy. However, that relationship doesn’t look clear in these data.
Could be a negative relationship between Orthographic and modality.au effect (on pv accuracy).
ppts %>%
pivot_longer(cols=c(Visual, Verbal, Orthographic, Manipulation), names_to="factor", values_to="score") %>%
ggplot(aes(x = score, y = wm.modalityau)) +
geom_point() +
geom_smooth(method="lm", formula="y~x") +
facet_wrap(facets=vars(factor)) +
theme_minimal() +
labs(
# title="IRQ vs PR Domain ∆",
# subtitle = "By-Participant",
# x = "IRQ Dimension Score",
# y = "PR Accuracy Domain ∆ (Visual - Social)"
)Simple linear models show no effect of either variable.
summary(lm(wm.modalityau ~ Verbal, data=ppts))
summary(lm(wm.modalityau ~ Orthographic, data=ppts))##
## Call:
## lm(formula = wm.modalityau ~ Verbal, data = ppts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6123 -0.9973 0.0761 1.1443 2.7429
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.05471 0.42143 0.130 0.898
## Verbal -0.03655 0.07276 -0.502 0.622
##
## Residual standard error: 1.885 on 18 degrees of freedom
## Multiple R-squared: 0.01382, Adjusted R-squared: -0.04097
## F-statistic: 0.2523 on 1 and 18 DF, p-value: 0.6216
##
##
## Call:
## lm(formula = wm.modalityau ~ Orthographic, data = ppts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3681 -0.6110 -0.0879 0.7227 2.8988
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.05471 0.39709 0.138 0.892
## Orthographic -0.13000 0.08128 -1.599 0.127
##
## Residual standard error: 1.776 on 18 degrees of freedom
## Multiple R-squared: 0.1244, Adjusted R-squared: 0.0758
## F-statistic: 2.558 on 1 and 18 DF, p-value: 0.1271
Exploratory.
Again, high orthographic loading correlates with a negative effect of pv-au (less accurate on auditory stimuli). Not sure why this would be the case
ppts %>%
pivot_longer(cols=c(Visual, Verbal, Orthographic, Manipulation), names_to="factor", values_to="score") %>%
ggplot(aes(x = score, y = pv.modalityau)) +
geom_point() +
geom_smooth(method="lm", formula="y~x") +
facet_wrap(facets=vars(factor)) +
theme_minimal() +
labs(
# title="IRQ vs PR Domain ∆",
# subtitle = "By-Participant",
# x = "IRQ Dimension Score",
# y = "PR Accuracy Domain ∆ (Visual - Social)"
)Simple linear models show no effect of either variable.
summary(lm(pv.modalityau ~ Orthographic, data=ppts))##
## Call:
## lm(formula = pv.modalityau ~ Orthographic, data = ppts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1403 -0.5573 -0.2812 0.4885 2.4288
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.14706 0.28736 0.512 0.615
## Orthographic -0.09539 0.05882 -1.622 0.122
##
## Residual standard error: 1.285 on 18 degrees of freedom
## Multiple R-squared: 0.1275, Adjusted R-squared: 0.07902
## F-statistic: 2.63 on 1 and 18 DF, p-value: 0.1222
Exploratory.
Orthographic shows the steepest relationship here too, but all seem to be pushed about by these two extreme points.
ppts %>%
pivot_longer(cols=c(Visual, Verbal, Orthographic, Manipulation), names_to="factor", values_to="score") %>%
ggplot(aes(x = score, y = `pv.modalityau:wm.modalityau`)) +
geom_point() +
geom_smooth(method="lm", formula="y~x") +
facet_wrap(facets=vars(factor)) +
theme_minimal() +
labs(
# title="IRQ vs PR Domain ∆",
# subtitle = "By-Participant",
# x = "IRQ Dimension Score",
# y = "PR Accuracy Domain ∆ (Visual - Social)"
)Simple linear models show no effect of either variable.
summary(lm(`pv.modalityau:wm.modalityau` ~ Manipulation, data=ppts))
summary(lm(`pv.modalityau:wm.modalityau` ~ Orthographic, data=ppts))##
## Call:
## lm(formula = `pv.modalityau:wm.modalityau` ~ Manipulation, data = ppts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8879 -0.3898 0.2314 0.7771 1.5262
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.15753 0.27271 -0.578 0.571
## Manipulation 0.03241 0.07578 0.428 0.674
##
## Residual standard error: 1.22 on 18 degrees of freedom
## Multiple R-squared: 0.01006, Adjusted R-squared: -0.04494
## F-statistic: 0.1829 on 1 and 18 DF, p-value: 0.674
##
##
## Call:
## lm(formula = `pv.modalityau:wm.modalityau` ~ Orthographic, data = ppts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7611 -0.2923 0.2208 0.7025 1.7898
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.15753 0.26437 -0.596 0.559
## Orthographic 0.06281 0.05411 1.161 0.261
##
## Residual standard error: 1.182 on 18 degrees of freedom
## Multiple R-squared: 0.06965, Adjusted R-squared: 0.01796
## F-statistic: 1.348 on 1 and 18 DF, p-value: 0.2609